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Abstract 



CO ' We consider the heavy-trafSc approximation to the GI /M / s queueing system in the Halfin- 

Whitt regime, where both the number of servers s and the arrival rate A grow large (taking the 
service rate as unity), with A = s — /3-\/s and /3 some constant. In this asymptotic regime, the 
, queue length process can be approximated by a diffusion process that behaves like a Brownian 

p ^ ' motion with drift above zero and like an Ornstein-Uhlenbeck process below zero. We analyze 

the first passage times of this hybrid diffusion process to levels in the state space that represent 
congested states in the original queueing system. 
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> 

^ ■ 1 Introduction 

o ■ 

■ Halfin and Whitt [8] introduced in their 1981 paper a new heavy-traffic limit theorem for the 

CN . GI/M/s system. They demonstrated how under certain conditions a sequence of normalized queue- 

length processes converges to a process that behaves like a Brownian motion with drift above zero 
and like an Ornstein-Uhlenbeck process below zero. We refer to this hybrid diffusion process as the 
Halfin-Whitt dijfusion. 

In [8] it is established that by setting the traffic intensity p = 1 — /?/ ^/s, (5 G (0, oo), the number 
of customers in the M/M/ s system can be roughly expressed as s -|- y/sX{t) for s sufficiently large 
. and [X{t))t>o the Halfin-Whitt diffusion. The boundary between the Brownian motion and the 

Ornstein-Uhlenbeck process can be thought of as the number of servers, and {X{t))t>Q will keep 
fluctuating between these two regions. The process mimics a single server queue above zero, and an 
infinite server queue below zero, for which Brownian motion and the Ornstein-Uhlenbeck process 
are indeed the respective heavy-traffic limits. As (3 increases, capacity grows and the Halfin-Whitt 
diffusion will spend more time below zero. 
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The diffusion process {X{t))t>o can tlius be employed to obtain simple approximations for the 
system behavior. The steady-state characteristics of the diffusion were studied in [8]. It is also of 
interest to study time-dependent characteristics like the mixing times, time-dependent distributions 
and first passage times to enhance our understanding of how the GI/M/s system behaves over 
various time and space scales. The mixing time is closely related to the spectral gap, which for the 
Halfin-Whitt diffusion {X{t))t>o has been identified by Gamarnik and Goldberg [6] building on the 
results of van Doorn [5] on the spectral gap of the M/M/ s queue. An alternative derivation of this 
spectral gap was presented in [12], along with expressions for the Laplace transform over time, and 
the large-time asymptotics for the time-dependent density. In this paper we derive results for first 
passage times to large levels. Such large levels typically correspond to highly congested states, in 
which users start receiving degraded service. An expression for the mean first passage time was 
derived in Maglaras and Zeevi [14]. We shall derive the Laplace transform of the first passage time 
density. From this Laplace transform, we can derive not only all moments, but also expressions for 
the first passage time density in various asymptotic regimes. 

Mathematically, determining the Laplace transform of the first passage or time-dependent dis- 
tributions for the present diffusion process involves analyzing a Schrodinger type equation with a 
piecewise parabolic potential function, or, equivalently, a Fokker-Planck equation with a piecewise 
linear drift. Such problems arise in a variety of other applications, such as linear systems driven 
by white noise [4, 3], the Kramers' problem [15] and escape over potential barriers [11]. Invariably, 
the solution involves the parabolic cylinder functions (see also [1, 13, 16] for more background on 
the parabolic cylinder function). The main results are presented in Section 2 and the proofs are 
given in Sections 3-5. 

2 Main results 

For the Halfin-Whitt diffusion process, define Tx{b) as the first passage time out of the interval 
(—00,6), starting at x < 6 with 6 > 0. Define 

(^,,6(e) = E[e-^^^(^)], (2.1) 

so that if P(x, t)dt = P (T^(6) G [t,t + dt]) 

POD 

^x,b{0)= e-^'Pix,t)dt, K(0)>O. (2.2) 
^0 

The first passage time density P satisfies the backward Kolmogorov equation 

Pt = A{x)Px + \B{x)Pxx; x<b, t>0 (2.3) 
with P{b, t) = 6{t) (the Dirac function) and 

We also require P and Px to be continuous at x = 0. Here the diffusion coefficient is B{x) = 1 + 
where c is the coefficient of variation for the interarrival distribution of the GI/M/s system. For 
GI=M we have c = 1, and in general we can rescale x so as to make B{x) = 2, which we henceforth 
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assume. 

Let Dy{z) denote the parabolic cylinder function with index v and argument z, which is defined, 
for example, by the integrals 



-z^/A roo 



r(-z^) Jo 

e 



i^/2TT Jc 



u 



(2.5) 
(2.6) 



Here, r(-) is the Gamma function, and the contour C in the second integral is a vertical Bromwich 
contour in the half-plane 3?(ti) > 0. It is well known that Dy{z) is an entire function of both index v 
and argument z, and various properties of Dy{z) are given in [1, Chapter 19] and [7, p. 1092-1095]. 
Define 



b 



M(M,5)^cosh(-7^^ i._.(-/3) 



(2.7) 



with D'_q{—/3) = —-^D_0{—P). Below we give expressions for (px,b{(^), where we must distinguish 
between the cases a; > and x < 0. 



Theorem 1. Let x <0. Then, with M{9;/3,b) as defined in (2.7), 

1 D_e{-li-x) 



fx,b{0) 



exp 



M(^;/3,6) D_e{-f^) 
Theorem 2. Let x > 0. Then, with M{9;(3,b) as defined in (2.7), 



/3(6 — x) 

2 ^ T 



(2i 



Using 



/3(x - b) 



sinh(fV/32 + 40) ^ sinh(^V-52 + 40 



sinh (|v//32 + 4^) M{6; 13, b) ^^^^ (^1^/32 + 45/ 



(2.C 



^[TAb)] = -^ip,,b{0) 



(2.10) 



we obtain after tedious calculations the following result for the mean first passage time, which is 
in agreement with the result obtained in a different manner by Maglaras and Zeevi [14]. 



Proposition 3. [14, Proposition 3] If x > then 



E[TAb)] = ^+(e^'-e^^ 



1 1 



If x < then 



E[T,(6)] = ^7^^ + 



3^«-«'/2dn. 



3^"-"'/2dn 



3u-v?/2 



1 



U 



du. 



(2.11) 



(2.12) 



/32 ' /5 

Though (2.11) and (2.12) are already fairly simple, we give some asymptotic formulas below 



3 



that yields further insight on the magnitude of the mean passage time (the derivation is standard 
and therefore omitted). 



Proposition 4. (a) For 6 — )• oo and b — x 

E[r,i 



oo, 



~ e 



1 1 



If b — X = 0(1) the above term should be multiplied by 1 
simplifies further to E [T^{b)] ~ V^^^-'^e^^e^^^'^ . 
(b) For P -oo and x,b = 0{\(3\) {possibly o(|/3|)), 



-t3{b-x) 



Ifb,(3 



(2.13) 
oo then (2.13) 



E[T,( 



b—x 

-13 ' 

^ + log ( 1 + § 



xG [0,6), 
X £ (— oo, 0]. 



(2.14) 



(c) For /3 oo with /3 = 0(6"^) and x = 0{b), 
E[T,(6)]- 



/3-2[e/^^-e'^^ + /3(x-6)], 2;/6g[0,1), 
/3-2[e^^- 1-/36], x/b<0. 



(2.15) 



We note that in (2.13) the mean first passage time is exponentially large and independent of the 
starting point x, in (2.14) it is asymptotically 0(1), and (2.15) represents the transition between 
these two cases, where E[Ta;(6)] = 0(6^). 

The Laplace transform ipx,b{(^) is analytic in the entire 0-plane, except for singularities in the 
range 5R(0) < 0. Hence, the asymptotic behavior of Tx{b) is determined by the singularity 0max 
closest to the imaginary axis. In fact, from Theorems 1 and 2 it follows that ^max will be the largest 
negative solution to 

i?_,(-/3)M(0;/3,6) =0. (2.16) 



Note that (2.7) and (2.9) are invariant under the change y^/?^ + 40 — )• — \/ /^^ + 40, so there is 
no branch point at = — /32/4. It seems impossible to find a closed-form solution to (2.16). We 
therefore consider several asymptotic regimes: 

(i) Large levels: 6 — )• oo and /3 fixed. 

(ii) Large levels and over/undercapacity: /3 — )• ±00, 6 — )• 00 at the same rate (|/3|/6 fixed). 

(iii) Small levels and undercapacity: /3 — )• —00 and 6 — )• 0. 

Regime (i) represents the situation of reaching highly congested states, corresponding to large 
levels 6. We have the following results: 

Proposition 5 (Regime (i)). If (3 < is fixed, and b — )• 00, then 



TT 

6^ 



l + ???4^ + 0(6-) 



bD' 



/32/4 



(-/?) 



(2.17) 
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If (3 = 0, and b — )• oo, then 



IT 

'462 



+ 0(6" 



(2.18) 



If 13 > is fixed, and b — )• oo, then 



1 + /3e/3V2 e-«V2dti 

J — oo 

We can generalize (2.18) to the case where 6 — oo and /5 — >• with 



1 



7 



+ ^(7) 



where w is the solution of the smallest absolute value to 

tan(-^/a;) tanh(-y/— w) 



-o; 



2 

7' 



(2.19) 

7 fixed, where we have 
(2.20) 

(2.21) 



It follows that if 7 = 0, a;(0) = 7r2/4 and then the leading term in (2.18) becomes a special case of 
(2.20). Also, if 7 = 2, uj{2) = 0, with w > for 7 < 2 and < for 7 > 2. 

From (2.19) we see that |0max| is exponentially small, which implies exponentially large time 
scales. While the result (2.19) is established analytically in Section 4, it can also be seen as a 
consequence of the following result. Let denote convergence in distribution. 

Proposition 6 (Exponential limit law). Let V be an exponential random variable with unit mean. 
Then, 



Ce-f^^T^ib) ^ V, as 6 ^00 



with 



00 \ — 1 







(2.22) 



(2.23) 



As mentioned in [14], Proposition 6 can be established using a limit theorem from the theory 
of regenerative processes. In Section 5 we give two proofs of Proposition 6. The first proof is 
probabilistic and uses the theory of regenerative processes as pointed out in [14], and the second 
proof is analytic and uses the exact expressions for the Laplace transform. 

We next give some results for the double limits 6 — >■ 00 with f3 — t- ±00, and also 6—7-0 with 
-00. 



Proposition 7 (Regime (ii)). If f3 



-00 and 6 — )• 00, then 
,2 



TT 



--0^ 

4^ 62 



2 



2 y/^ Ai(o) 

^) aF(o) 



(2.24) 
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where Ai(-) is the Airy function, /f /3 — )• oo and b — )• +00, then 



/3 



Proposition 8 (Regime (iii)). )• with (3 

1 



4 



e-'-^e-^'/^^ (2.25) 
—00, we let Bt, = b{—l3)^^^ and for fixed, 

(2.26) 



2) ^' 



where rj = r}{B^) is the minimal solution to 



-T] cot 



i?*2-l/^ 



Ai{ri) • 



(2.27) 



We comment that (2.24) remains valid for /3 — )• —00 with 6 > fixed, (2.25) remains vahd for 
/? — 7- 00 with 6 > fixed and can be obtained as a hmiting case of (2.19) (for (3 — t- 00). 

When 5* 00 it follows from (2.27) that rj ^ with r] ~ -2'^/^it^B:;^, while if B^ 0+ we 
have r] ro = —2.338 . . . where rg is the least negative root of Ai(z) = 0. If r/ = < where is 
the least negative root of Ai'(z) = 0, then B^ = i7r2^/^(— r*)~^/^. 



3 Proofs of Theorem 1 and Theorem 2 

From (2.3) and (2.2) it follows that the Laplace transform (px,b{G) = Q{x;9) satisfies the ODE 

Qxx + A{x)Qx = eQ, x<b, (3.1) 

with the boundary condition Q{b;9) = 1 and the interface conditions Q{0^;9) = Q{0^;9) and 
Qx{0^', 9) = Qx{^ ] 9). For X > we have A{x) = —(3 and then (3.1) admits solutions of the form 
e"^' where a = a± = ^[/3ib /J^ + 40]. For x < 0, ^4(3;) = —x—f3 and then (3.1) becomes the Hermite 
equation, and the only solution that decays as x — )• —00 is proportional to e^^^^e^^^'^D^g{—x — j3). 
It follows that 

Q(x;g) = fcoexp(^- + ^J % ^ < 0, (3.2) 

and 

Q(x; 9) = A;ie°+(^)(^-^) + (1 - A;i)e°-(^)(^-''), < x < 6, (3.3) 

where and ki are independent of x. Here we wrote the solution in (3.3) in such a way so to 
automatically satisfy the boundary condition Q{b\9) = 1. To determine /cq and ki we can use the 
interface conditions at x = 0, which imply that 

feo = feie-"+(^)'' + (1 - fci)e-"-(^)^ (3.4) 
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and 



D' 



1 



D_e{-P) 2 



+ AO [fcie-"+(^)* - (1 - A;i)e-"-(^)^ 



(3.5) 



Solving the algebraic system in (3.4) and (3.5) for /cq and ki leads to the expressions in Theorems 
1 and 2. 



4 Brief derivation of Propositions 5, 7 and 8 

We discuss the various asymptotic formulas for ^max- The expressions follow from routine manip- 
ulations of the parabolic cylinder functions that appear in (2.7) and (2.16). Consider the solution 
to D_Q{—j3)M{9; /3,b) = which is equivalent to 

D_0{-/3)cosh Qy^^TI^) = -^=^=Z)'_,(-/3)sinh Q^^^TI^^ . (4.1) 

We analyze this transcendental equation in various limits and find the least negative root ^max) 
only sketching the main points in the calculations. 

For 6—7-00 and /3 > we use the Taylor expansion, for 0—7-0, 



D 



=!!_ Qji^n 
-P) 2 



e-"'/2cln + 0(^2) 



and, for — )• and h — )- 00, 



coth(^V/52 + 40) = l + 2e-^^[l + O(0,e-''^)]. 



Also, + AO = (2//3)[l - 20//32 + 0{e'^)], so that (4.1) is equivalent to 

./3 



and hence 



72 



J —00 



l + 2e-'"' + O(0,e 



-26/3 \ 



(4.2) 



(4.3) 



(4.4) 



-29 



732 + 



2e" 



(4.5) 



which leads to the exponentially small ^max in (2.19). A similar calculation applies for fixed 6 > 
and (3 — ?- 00, which leads to (2.25) in Proposition 7. 

For 6—7-00 with /3 < the solution ^max will be close to the apparent branch point at = — /32/4. 
From (4.1) if we set 9 = - fi"^ / A: - ^/b'^ we obtain 



+ 0{9 + /32/4) 



tan(\/j7). 



(4.6) 



D-e{-P) — ' ' ' 

If 6 —7- 00 with a fixed /? < 0, V^l must be close to a zero of the tangent function, so that VL ~ 7r2/4 
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(for the least negative solution ^max)- Then estimating the difference — 7r/2 using (4.6) leads 
to (2.17). 

If 6 — )• oo with /3 — 7- we again obtain (4.6) and for small (3 we can further approximate 

D'_e{-P) ^ 2 ^ 26 

Then (2.20) follows from (4.6) with VL replaced by cj = b^O^g.^ ~ 7^/4- If /3 = we can express the 
parabolic cylinder functions D_g{0) and D'_g{0) in terms of Gamma functions, and (4.1) becomes 

tanh(67g) ^ -r(|) 

^/2r(^)' ^ • ^ 

Then (2.18) follows by solving (4.8) for 6'niax with b ^ oo, where 6'max n'^/iib'^) = 0{b~'^). 

Then r(^) ~ and r(|) ~ |. 

If /3 — 7- — oo and 6 — )■ oo or 6 = o(l), (2.24) follows by a calculation similar to (4.6), except 
now Dp2/4^{—f3) becomes proportional to Ai(0) in this limit. Finally, to obtain (2.26) we use the 
approximation 

-^D'_,{-(3) f-^^y/'Al'irj) g + /3V4 ^^^^ 



D_e{-f3) V 2 ; Ai(r/) ' ' (-/3/2)2/3 

which applies for -(3 ^ oo and -6^ = /374 + 0(|/3|2/3). With (4.9), (2.27) follows from (4.1). This 
concludes the sketched derivation of Propositions 5, 7 and 8. 



5 Two proofs of Proposition 6 

5.1 Probabilistic proof 

The crucial observation is that the one-dimensional Halfin-Whitt diffusion process is an ergodic one- 
dimensional diffusion process that has, by the strong Markov property, the origin as a regeneration 
point. To make this formal, let y be a fixed positive number, where y < b. This value can be used 
to construct a sequence of random times {Tk}k>i, where 

Ti = mi{t > : X{t) = 0, sup X{s) > y} (5.1) 

0<s<t 

and for each n > 1, 

T„+i = inf{t > r„ : X{t) = 0, sup X{s) > y}. (5.2) 

Tn<S<t 

The Halfin-Whitt diffusion is a regenerative process with respect to these regeneration times, which 
form a delayed renewal process that, with probability one, has a finite number of points in each 
compact interval. To make matters simpler to state, we assume that X{0) = 0, but the procedure 
outlined here can be adjusted for any arbitrary initial condition. 

Define a{b) to be the probability that our process gets above level b in the random interval 
[0,Ti]. Due to the Strong Markov property, along with the fact that y < b, we observe that this 
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probability is just the probability that a Brownian motion with drift — /3 and diffusion coefficient 2 
reaches level b before level 0. Hence, 



aib) 



1 



(5.3) 



Moreover, the expected length of each regenerative cycle is just the expected amount of time it 
takes the Halfin-Whitt diffusion to reach level y, starting from 0, plus the expected amount of time 
it takes the diffusion to go from level y back to 0. In other words, 



E [Ti] = {e 



I3y 



I) 



1 1 



Therefore, by [2, Theorem 4.2 on p. 181], we see that as 6 — t- oo, 

a{b) 



T.{b) V 



(5.4) 



(5.5) 



where V is an exponential random variable with rate one. Moreover, observe that for each y < b 

a{b) 1 



E[Ti 



lib 



1) 



(5.6) 



which does not depend on y. Hence, as 6 — )■ c«, we obtain (2.22). A very similar asymptotic result 
carries through as well for E [Tj(5)], due to Asmussen [2, Proposition 4.1 on p. 180]. From (2.22) it 
follows that the time it takes to reach a level b is roughly exponential in 6, which says that extreme 
congestion is not observed on relatively short time scales. For /3 > the mean first passage time is 
E \Tx{b)] ~ l/|^max| for X bounded away from b and 6 — oo. 



5.2 Analytic proof 

Consider the Laplace transforms in (2.8) and (2.9) on the scale 9 = O(|0max|) = 0(e~ 
and then scale time as t = T/\9raax.\ = 0(e*^). Since L'o(— /3) = e~'^ we have 



exp 



Then we write M in (2.7) as 



M 



sinh 1^/32 + 40 



V/S^ + 49 



VP^+A9 coth i^- + A9 



D-e{-p) 



in (2.19), 
(5.7) 



(5.8) 



For —7- we use (4.3), and the recurrence relation for the parabolic cylinder function Dp{z) 
— ^Dp{z) + pDp^i{z), which for p — )• yields 



Dp{z) 



Z D_l(z) ^. n. 



(5.^ 
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Then (5.8) becomes 



M = —J'P/^ 



2/3 
1 



b/3/2 



-b/3/2 



(1 + 2e-^^ + 0(06" 

~l~ l^maxl 



+ e'^'/4D_i(-/3y 



+ 2/3e 



))(/? + ^ + 0{e^)) - /3 + 20e'^'/^I)_i(-/3) + O(02) 
-bp 



(5.10) 



^ = Od^inaxl), as -D-i(-/3) = e/^'/4 /f^ e~"'/2du = e'^'/^ e^^e-^'/^du. Then inverting the 
transform in (2.8), using (5.7) and (5.10), yields (with Br a vertical Bromwich contour in the 
half-plane 3ft(^) > 0) 



2TTi 



Br 



1+C 



(5.11) 



since 9t = \i 9 = l^maxlC- This yields the exponential limit law for x < 0. 

If X > so that both 6 — )• oo and h — x ^ oo, then an analogous calculation using (2.9) again 
leads to (5.11). However, if x, 6 — )• oo in such a way that b — x = 0(1), then (2.9) leads to the limit 



and this inverts to, on time scales t 

p(x,t)~l! 



+ 



+ 



1 



,{x-b)(3 



(5.12) 



0(e' 



b|3^ 



27rz 



,(x-fe)/3 

Br 

{x-b)l3^^rp^ 



.{x-b)P 



+ 



-e( 



(5.13) 



The term proportional to 5{T) represents a probability mass on the large time scale, which cor- 
responds to sample paths that hit 6 in a short time (at least t = o(e^^)). The actual density of 
P{x,t) does not have mass at t = 0, and finding the expansion of the density on shorter time ranges 
would require a different asymptotic analysis, and a different approximation to M. We refer the 
reader to [9, 10], where such problems were analyzed in detail for some simpler models. There both 
exponentially large and t = 0(1) time scales were considered, and they were related to one another 
by asymptotic matching. 
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